Theoretical Population Biology

A sudden reduction in population size increases the rate of genetic drift, reducing variability and increasing the mean level of homozygosity. The resulting increased exposure of recessive or partially recessive, strongly deleterious alleles to selection against homozygotes may lead to their being purged from the population, potentially allowing mean fitness to increase after an initial decline, and accelerating the decline in inbreeding depression associated with reduced variability. However, detailed population genetic theory on the effects of population bottlenecks on mean fitness and inbreeding depression remains limited. We develop a theoretical framework for small, randomly mating populations founded from a large population near mutation-selection-drift equilibrium, using both simulations and approximate analytical predictions. These provide quantitative predictions for the dynamics of the populations mean fitness and level of inbreeding depression following a bottleneck. In particular, we derive an approximate expression for the time needed for mean fitness to recover after an initial decline; such a recovery requires selection to be sufficiently strong relative to drift and mutations to be sufficiently recessive. In contrast, weakly deleterious mutations cause reductions in mean fitness and inbreeding depression that are similar in size to those predicted from increases in neutral homozygosity.

Darwinian fitness is equated here with invasion fitness and defined as the quantity determining the fate--certain extinction or possible spread--of a single mutant type. We derive it, together with its phenotypic derivative, for evolution in group-structured populations under limited genetic mixing, where the demography of the focal species and its environment is modeled as a discrete-time stochastic process. Reproduction, physiological development, dispersal, and survival are influenced by interactions within and between groups and by environmental fluctuations within and across generations. Using multitype branching processes in random environments, we show that invasion fitness is predicted by a stochastic growth rate that can be represented biologically in two meaningful genealogical ways. First, as the long-term geometric mean of the expected per-capita number of mutant copies produced per time step by a representative member of the mutant lineage. Second, as the the long-term geometric mean of the expected reproductive-value-weighted per-capita number of mutant copies produced by such an individual. This latter representation is useful for computing the phenotypic directional derivative of invasion fitness. Moreover, this derivative can be written as an actor-centered inclusive-fitness effect derived from properties of the resident population process. This effect depends on class-specific fitness differentials, relatedness, reproductive values, and class frequencies. However, unless generation- and class-specific fitness defines a stochastic matrix, the derivative does not separate stochastic reproductive values from relatedness and class frequencies, and must be evaluated by simulations. In summary, we formalize invasion fitness biologically quite generally and show how Hamiltons marginal rule is deduced from it.

We study a polygenic trait under stabilizing selection at statistical equilibrium, where genetic effect, mutation rate and mutational bias are heterogeneous across loci. The model assumes L biallelic sites subject to reversible mutations, each allele described by its frequency in the population. Using a diffusion approximation, a mean-field approximation and neglecting linkage disequilibrium, we predict consistent phenomena across several regimes of selection: (1) a small deviation {Delta}* of the trait mean from its optimal value appears and persists due to genetic mutations not aligned with selection; (2) while this deviation is often undetectable at the trait level, it leaves a substantial signature at the locus level by favoring alleles reducing it, resulting in genic selection with mean coefficient s* proportional to -{Delta}* acting pervasively; (3) with stronger selection on the trait, (3a) the value of {Delta}* is decreased but the intensity of genic selection is increased in inverse proportion, resulting in an essentially constant, non negligible value of s*. We show how the stationary distribution of allelic frequencies can be obtained from {Delta}*. The latter can then be characterized as the solution to a fixed-point equation. Finally, we quantify several macroscopic observables of interest (genetic variance, description of the fluctuations of the trait mean as an Ornstein-Uhlenbeck process). The orders of magnitude of the macroscopic observables can be derived on a wide region of the parameter space. The model shows good fit and can straightforwardly be extended to accommodate pleiotropy, dominance, and some forms of epistasis. We also discuss the different breakdown which may occur (Bulmer effect, Hill-Robertson effect, breakdown of the Ornstein-Uhlenbeck approximation for the dynamics of the trait mean, depletion of genetic variability due to low mutation rates).

Price equation and genotype dynamics are two methods for studying the fixation of one allele by natural selection in a diploid population. There are two strict monotonicity conditions that imply the fixation of one allele. The genotype dynamics is called Haldane monotone if the relative frequency of one allele strictly increases along all solutions of the genotype dynamics, so this allele is fixed. In this paper, we show that the genotype dynamics is Haldane monotone if and only if the right-hand side of the Price equation is always strictly positive. The other strict monotonicity condition requires that the relative frequency of a homozygote strictly increase according to the genotype dynamics. For example, in a model where the genotype dynamics is governed by interactions between individuals, the cost-accepting homozygote is fixed by natural selection if the other genotypes always receive a smaller average gain from all interactions than the cost-accepting homozygote. Both monotonicity conditions require that the interaction is not well-mixed in the population. These two conditions are not equivalent. In addition, we give a non-monotonicity condition, which also implies the fixation of a homozygote. The fixation of a homozygote depends on the phenotypic payoff of the interaction, the genotype-phenotype mapping, and the interaction scheme. In a sexual population, the interaction scheme of siblings depends on the mating system, and so do the conditions of fixation of the cost-accepting homozygote. We present examples showing that if we only change the monogamous mating system, assuming panmixing or mating assortativity, then the condition for the fixation of the cooperator homozygote is b > 2c and b > c, respectively.

Interest in quantifying linkage disequilibrium (LD, non-random associations of alleles at different loci) has skyrocketed in recent years as researchers have focused on use of LD in genome-wide association studies (GWAS), for studying historical demography, and for estimating effective population size (Ne). The most widely used LD metric is r2 = the squared correlation of alleles at a pair of loci. Despite a half century of efforts, developing an unbiased expectation of r2 as a function of the many factors that can affect it (physical linkage, genetic drift, selection, migration, mutation, mating systems) remains elusive. Furthermore, even when all of these other factors are absent, empirical estimates of r2 are upwardly biased by sampling a finite number (S) of individuals, and that must be accounted for if one wants to focus on the desired signal of LD. Previous approaches to estimate [Formula] have been shown to be biased to greater or lesser degrees. The purpose of this short paper is to demonstrate that a simple and apparently exact expression for [Formula] does exist for the special case where sampling error is the only factor contributing to r2, in which case [Formula] = 1/(S - 1). When other factors contribute heavily to LD, [Formula] shrinks toward 0 as empirical r2 [->] 1. However, for estimating contemporary Ne with unlinked markers, empirical r2 will generally be small and 1/(S - 1) will provide a robust estimate of [Formula].

Understanding how mutations evolve on Y chromosomes is central to explaining the origin, diversity and persistence of sex chromosomes. Mutations occurring on the Y chromosome in sexual populations experience selective dynamics that differ markedly from those on autosomes, due to a reduced effective population size and the presence of large non-recombining regions containing alleles maintained in a permanently heterozygous state. These specific features alter gene transmission in the Y chromosome population compared to autosomes, even within the same pedigree. Here, we provide a two-sex diploid Wright-Fisher model that explicitly incorporates both sex chromosomes and autosomes within a unified population framework, in order to capture the influence of these specificities on the fate of mutations, not only considering fixation probabilities but also segregation times. We use diffusion approximations and provide analytical and numerical tools to compute these quantities across a wide range of parameters and selection regimes. We recover classical results on fixation probabilities in various scenarios, including purely beneficial, deleterious or overdominant mutations, and extend them in the light of mean segregation time, a key but often overlooked determinant of evolutionary outcomes over finite timescales. In particular, our analyses show that overdominant mutations are overall more likely to fix in observable time windows on the Y chromosome than on autosomes. Individual-based simulations corroborate our approximations and highlight parameter regimes where the theoretical approach is particularly useful, especially for parameter values inducing long segregation times or small fixation probabilities, for which simulations are impractical. Our results provide a comprehensive and tractable framework for clarifying how chromosome-specific features shape evolutionary dynamics beyond fixation probabilities alone.

We study one of the simplest scenarios of polygenic selection that can be imagined: a subdivided population of diploid individuals expressing an additive trait under spatially homogeneous stabilizing selection. We are interested in the amounts of variation that can be maintained at mutation-selection-migration-drift equilibrium, at individual loci and at the level of the trait, within and among subpopulations. We derive analytical approximations for variance components and summary statistics such as FST and QST under the assumptions of the infinite-island model and compare these with individual-based simulations. We find that: (i) There is a critical migration threshold (which depends on effect sizes of trait loci) below which population structure strongly inflates genic variance in the subdivided population to levels well above those in a panmictic population. Variation within each subpopulation is maximized close to the critical migration rate. (ii) The genetic basis of trait variation across subpopulations is most similar close to this migration threshold and (counter-intuitively) decreases for higher migration rates. This has consequences for the portability of Genome-Wide Association Studies (GWAS) between subpopulations, i.e, the extent to which loci with large contributions to variance in one subpopulation explain variance in other subpopulations. (iii) An analytical mean-field approach based on the single-locus diffusion approximation, together with effective migration and selection parameters (to account for associations between loci), very accurately predicts various quantities.

Minimum viable populations (MVPs) are population levels large enough to surmount risk from demographic, environmental, and genetic stochasticity. MVPs are estimated by biologists to guide conservation practices. However, MVPs are generally estimated for a target population without regard for how they interact with intra- and inter-species population dynamics in the broader ecological community. Thus, how and why population dynamics interact with MVPs imposed by conservation biologists remain unclear. When MVPs are imposed on a continuous population model, traditional analyses fail to capture the range of possible outcomes those MVPs create. Here, we describe viability space decomposition (VSD) as a mathematical tool to systematically analyze the potential crossing of MVPs during population dynamics. We demonstrate that different extinction and survival outcomes can be recovered from a model with imposed MVPs using three VSD concepts in junction with a traditional phase portrait: mortality manifolds which separate conditions that lead to different existential outcomes, ordering manifolds which determine the order of extinction events for multiple populations, and collapse manifolds which determine the survival or extinction of one species given the loss of another. We employ these methods with a standard consumer-resource model, and the methods can be scaled to systems with more species. VSD is a useful tool for conservation biologists and community ecologists concerned with boundary crossing problems in any dynamical system.

1Although resources are typically distributed continuously in space, species distributions often organize into discrete clusters. In his seminal paper [36], Turing demonstrated that such clusters can spontaneously arise in population densities, even when populations evolve in environments with continuously varying conditions. This phenomenon is known as Turing instability. In this work, we focus on two models grounded in population dynamics: a one-dimensional model based on the nonlocal Fisher-KPP equation, and a two-dimensional model involving an environmental gradient. We show that phenotypic clusters (sometimes referred to as "species") emerge in these models. We prove that they do not emerge because of Turing instability, but because of stochasticity, and that they disappear when stochasticity is reduced. First, for both models, we start our simulations with initial populations uniformly distributed in the state space. We show that phenotypic clusters quickly emerge and that the distances between them depend on the population size, that is, on the degree of stochasticity. Next, we start from already clearly defined phenotypic clusters. We identify three regimes in the connection between population size, the initial distances between clusters, and the distances between clusters at equilibrium. Last, on the two-dimensional model, we relax the hypothesis of complete clonality by varying the effective recombination rate, explore its effect on phenotypic clustering, and show that phenotypic clustering decays drastically with slight recombination.

Alleles with opposing effects on fitness characters are said to exhibit selectional antagonistic pleiotropy (broadly construed so that effects are not necessarily confined to the same individual). A number of theoretical investigations considered the case where a pair of alleles at a locus influences two fitness components and derived the conditions giving rise to stable polymorphism under various assumptions about the mode of trait-interaction. Strikingly, many of these analyses concluded that the potential for maintaining polymorphism is strongly constrained by the joint influence of two factors: (1) the prevalence of weak selection coefficients over coefficients of large magnitude, and (2) the absence of beneficial dominance reversals (where the deleterious effects of each allele are partially or completely masked in the heterozygous genotype). Consequently, the conclusion that selective polymorphism is unlikely to be maintained by intralocus mechanisms of antagonistic pleiotropy has achieved widespread acceptance. Here we argue that such conclusions do not apply to any of the following models of antagonism: (i) additive trait-interaction, (ii) multiplicative trait-interaction, (iii) bivoltine selection, (iv) soft selection, (v) hard selection, and (vi) sexual antagonism. We demonstrate that the parameter space giving rise to stable allelic variation is quite large throughout, and moreover, the plenitude of suitable parameters neither depends on the strength of selection nor requires dominance reversal. Dominance coefficients associated with stringent conditions for stable polymorphism are shown to be atypical as compared to all feasible parameters, and best regarded as an outcome of adherence to a special relation: dominance with a constant magnitude and direction, which includes the case of additive allelic effects at a locus. Properties of single-locus equilibria (heterozygosity, allele frequency differentiation) are investigated, as well as the contribution of dominance schemes to the genetic variance in fitness characters in populations at multilocus linkage equilibrium. Author summaryAllelic variants at a locus with opposing effects on multiple fitness components (antagonistic fitness pleiotropy) have long been appreciated as a possible source of balancing selection. The prevalence of polymorphism owing to this form of natural selection, however, has been doubted on theoretical grounds due to the fact that standard assumptions of genetic models (namely, constant magnitudes for the dominance coefficients) are hardly conducive to the maintenance of polymorphism. The major exception to this conclusion lies with schemes that exhibit dominance reversal (where the direction of dominance for antagonistic alleles flips across fitness components). Here we conduct a geometric analysis of the space of polymorphism-promoting dominance parameters and conclude that the conditions for maintaining balanced alleles is unrestrictive, with non-reversals playing an underappreciated role.

Mark-recapture analyses on the delineation of natural populations between areas often assume random sampling, with a between/within (B/W) area resighting ratio that declines towards zero as the population components of two areas become more-and-more isolated from one another, with fewer-and-fewer individuals mixing between areas. I use an individual based population model split in two areas to simulate this result, analysing also for the potential effects of the space-time fidelity of the mark-recapture sampling in the areas. I find that small B/W resighting ratios--that traditionally is taken as evidence of population isolation--can easily be observed within a completely mixing population if a random sampling scheme is restricted in space and/or time. Random sampling within restricted areas and time windows is not sufficient to estimate mixing rates and population isolation between areas, unless the resighting rates are analysed by a method that accounts both for the space-time fidelity of the scientific sampling scheme and the space-time fidelity of the distributional behaviour of the individuals in the population.

Most eukaryotes reproduce using partial clonality, for which appropriate population genetic models remain limited. This gap constrains our ability to accurately reconstruct past population dynamics, predict future trajectories, and infer the evolutionary processes involved. We present a Wright-Fisher-like model tailored for tracking the mean and the variance of genotype frequencies over generations at one locus with multiple alleles in a same finite-sized population with mutation. Different initial conditions and rates of clonality generate unique mean trajectories of genotype frequencies. Partially clonal populations converge to the same unique stable equilibrium as exclusively sexual populations, that only depends on the reciprocal mutation rates between alleles. The dynamics unfold in two phases: First, genotype frequencies move towards Hardy-Weinberg proportions; Then iterate along the Hardy-Weinberg proportions until reaching the stable equilibrium. Mean allele frequencies and gene diversity remain unchanged by different rates of clonality along the trajectories. Instead, clonality influences the speed at which populations return to Hardy-Weinberg proportions and thus shapes the temporal sequence of genotype frequency distributions over generations. Variance around each mean trajectory depends only on parental genotype frequency distributions and population size, not on clonality. Taken together, these explain why both negative and positive Fis values are expected in partially clonal populations, and why variance of Fis across loci is a reliable proxy for inferring clonal rates. Our model will enable the analysis and prediction of changes in genotype frequencies within monitored populations, and will support future inference methods relying on time-series genotyping data from a target population. HighlightsO_LIOut of equilibrium, sexual and clonal populations share the same two-step dynamics. C_LIO_LIFirst, return to Hardy-Weinberg parabola impacted by rates of clonality; Then, iteration along this parabola until reaching equilibrium that only depends on mutation rates C_LIO_LIIncreasing clonality change the speed and direction of mean dynamics out of Hardy-Weinberg parabola without affecting mean allele frequencies C_LIO_LIVariance around mean dynamics depends on parental genotype frequencies and population size but not affected by clonality C_LI Graphical abstract O_FIG O_LINKSMALLFIG WIDTH=200 HEIGHT=98 SRC="FIGDIR/small/717696v1_ufig1.gif" ALT="Figure 1"> View larger version (13K): org.highwire.dtl.DTLVardef@1207e9dorg.highwire.dtl.DTLVardef@587d2dorg.highwire.dtl.DTLVardef@18224eborg.highwire.dtl.DTLVardef@145e2ed_HPS_FORMAT_FIGEXP M_FIG C_FIG

We introduce a spectral existence criterion for the evolution of cooperation in the form of the inequality{lambda} maxb > c, where{lambda} max is the leading eigenvalue of an interaction operator encoding population structure, and b and c represent benefit and cost tradeoffs, respectively. Nowaks five rules for the evolution of cooperation correspond to cases in which the cooperation condition reduces to a scalar assortment coefficient. These results follow from the Price equation, which sheds light on a long-standing debate on the role of inclusive fitness and evolutionary dynamics in explaining the evolution of cooperation.

We extend previous models of cumulative cultural evolution by incorporating structured populations and social networks. We examine how connectivity and network topology shape the accumulation of cultural complexity under unbiased (copy randomly), indirectly biased (copy successful individuals), and directly biased (copy successful traits) transmission. We consider random, scalefree, and small-world networks, as well as the communication structures introduced by Mason and Watts, and derive analytical approximations for the homogeneous case. We find that the effects of social structure depend strongly on the form of transmission bias. Under unbiased transmission, network effects are weak except at very low connectivity. Under indirect bias, cultural complexity increases with connectivity, whereas direct bias shows optimal performance at intermediate connectivity, reflecting a trade-off between diffusion and diversity. Differences across topologies are generally modest once the average degree is fixed. Overall, our results show that no single social structure universally promotes cumulative cultural evolution; instead, its effects depend primarily on the dynamics of learning and innovation.

The basic reproduction number R0 is central to malaria epidemiology, yet it is typically treated as a static quantity derived under memoryless assumptions for mosquito demography. In natural systems, however, mosquito populations are shaped by delayed processes such as larval development and density-dependent feedback, introducing biological memory into vector dynamics. We develop a minimal delay-based framework that incorporates this memory into the Ross-Macdonald model by describing adult mosquito abundance with a retarded differential equation. This formulation induces a time-dependent transmission potential R0(t). Using complex analysis and the argument principle, we derive an explicit stability threshold [Formula], which separates stable from oscillatory transmission regimes. Near this threshold, delayed feedback produces slow relaxation times and sustained transient oscillations, implying that transmission potential may vary intrinsically even in the absence of external forcing. To account for ecological variability, we extend this deterministic condition into a probabilistic framework and define the stability probability as [Formula]. Numerical simulations and global sensitivity analysis show that recruitment and developmental delays are the primary drivers of instability, while adult mortality has a weaker stabilizing effect. These results indicate that malaria interventions may influence not only the magnitude of malaria transmission but also its dynamical stability. By linking delay dynamics, transmission theory, and uncertainty quantification, this framework provides a basis for stability-aware modeling and interpretation of malaria transmission under ecological variability. Author summaryMalaria transmission is often summarized by a single number, R0, treated as a fixed indicator of whether transmission will increase or decline. This assumes mosquito populations respond instantly to environmental conditions. In reality, mosquitoes develop through stages where larval conditions, such as crowding, nutrition, or temperature, affect adult populations only after a delay. This creates biological memory: todays mosquitoes reflect past environments. We show that this memory can fundamentally reshape transmission dynamics. When developmental delays are included, transmission potential is no longer constant but can fluctuate over time, even in stable environments. These fluctuations can persist or amplify depending on the balance between mosquito growth, mortality, and delay. As a result, variability in mosquito abundance or malaria transmission may arise from intrinsic dynamics rather than external drivers alone. Under ecological variability, stability becomes probabilistic, allowing estimation of how likely transmission is to remain stable. Interventions that reduce larval productivity or increase adult mortality may therefore both lower transmission and make it more predictable, improving interpretation and control strategies.

Melanoma is a cancer of the melanocyte, known to have an ability to readily switch between different transcriptional cell states that convey different phenotypic properties (e.g. hyper-differentiated, neural crest-like). This ability is believed to underpin intratumour heterogeneity and plastic adaptation, which contributes to resistance to therapy and immune evasion of the tumour. Therefore, understanding the mechanisms underlying acquisition of transcriptional cell states and cell-state switching is crucial for the development of therapies. We model a minimal gene regulatory network comprising three key transcription factors, whose varying gene expression encodes different melanoma cell states, and use deterministic spatiotemporal differential-equation models to study gene-expression dynamics. We exploit an approximation, based on cooperative binding of transcription factors, in which the models are piecewise-linear. We classify stable states of the local model in a biologically relevant manner and, using a naive model of intercellular communication, we explore how a population of cells can take on a shared characteristic through travelling waves of gene expression. We derive a condition determining which characteristic will become dominant, under sufficiently strong cell-cell signalling, which creates a partition of parameter space.

Gene drives, genetic constructs that can spread deleterious alleles in wild populations, have the potential to address some of the major pressing challenges of the Anthropocene such as invasive species, spread of disease vectors, and agricultural pests. However, responsible and effective deployment of gene drive requires taking into account the complex nature of real-world population connectivity networks. In particular, it is unclear how the topological position of the deployment site affects the spread process and its final outcome. Here we develop a framework for modeling gene drive spread in population connectivity networks, and study the eco-evolutionary dynamics of gene drive spread under complex population structures. We investigated the relationship between the position of the deployment site in the topology of the network and whether the gene drive is eventually lost, fixed, or maintained at an intermediate frequency. We identified network centrality measures of deployment sites that are highly correlated with the outcome of deployment for different gene drive designs and across diverse network topologies. We also show that there is a trade-off between the time-to-fixation and the final outcome, implying that multiple centrality measures of the deployment site would need to be considered when aiming to achieve rapid and successful population control using gene drives.

The practical utility of many modern phylogenetic comparative methods can depend on how accurately mathematical models capture the evolutionary process of traits. Boucher and Demery (2016) described a new quantitative trait model, Brownian motion with reflective limits, that they anticipated might be of use in testing hypotheses about a particular sort of constraint on phenotypic character evolution. Since their analytic solution for the probability function under this bounded evolutionary scenario was not practical to evaluate for reasonably-sized trees, Boucher and Demery (2016) also identified a creative technique for computing the likelihood of their model. The basis of this methodology derives from the convergence of an equal-rates, symmetric, ordered Markov chain and continuous stochastic diffusion in the limit as the number of steps in our chain goes to {infty} (or, alternatively, as their widths decrease towards zero). We refer to this convergence in the limit as the discretized diffusion approximation or (more compactly) the discrete approximation. We realized that this discrete approximation of Boucher and Demery (2016) unlocked a number of additional models for the phylogenetic comparative analysis of discrete and continuous trait data, and we explore several of these in the present article. Specifically, we examine application of this discretized diffusion approximation to the threshold model from evolutionary quantitative genetics, to a new "semi-threshold" trait evolution model, to a joint model of discrete and continuous traits in which the discrete trait influences the rate of evolution of our continuous character, as well as a model where precisely the converse is true, and to a discrete character dependent multi-trend trended continuous trait evolution model. We conclude with some context for the origins of our article and discussion of other possible applications of this powerful approach.

Article summaryGene interactions are common, yet additive genetic models often predict short-term evolution and breeding response. This study argues that additivity can arise because populations sample only a small neighbourhood of a curved fitness landscape. In additive channels, genetic variation is small enough that local curvature contributes little to heritable fitness differences. The study defines an additivity index ([A]g) that compares variance from the local slope of log-fitness with variance from curvature, and links this ratio to expected prediction accuracy under Gaussian assumptions. A selection-inheritance framework shows when additive channels persist and when populations leave them. It yields testable predictions.

Animals such as bees, ants, wasps, termites, and naked mole-rats live in colonies in which a single queen is the only female reproductive, an arrangement known as eusociality. Eusocial animals are known for their remarkably long lifespans. It has been argued that longevity becomes selected when queens are shielded from "external mortality". While such protection may contribute, we find a deeper reason: the eusocial reproduction strategy itself inherently creates selection for long lifespans. Lifespans typically reflect two processes: the baseline risk of death and the rate at which this risk increases with age. Each is a parameter in the Gompertz mortality equation. We show that the mathematical properties of eusocial reproduction lead to slowly-growing, older populations where selection acts more strongly on the rate at which risk increases than on the baseline risk. In addition, we show that channeling reproduction through a single female also selects for longevity, which we term the "queen effect". Thus, the dynamics of eusocial reproduction select for longer lifespan. More broadly, these results show that reproductive structure and population growth dynamics can fundamentally shape selection on lifespan, with implications outside eusocial systems as well.